This code reproduces the trial-level analyses reported in the following manuscript:

Cosme, D., Mobasser, A., & J. H. Pfeifer. If you’re happy and you know it: Neural correlates of self-evaluated psychological health and well-being

load packages

if(!require('pacman')) {
    install.packages('pacman')
}

pacman::p_load(tidyverse, gtools, GGally, sjstats, lme4, lmerTest, knitr, ggeffects, kableExtra, psych, install = TRUE)

define aesthetics

rois = c("#006989", "#56445D", "#8EC922")
constructs = c("#FEC601", "#254E70", "#F43C13")
instructions = wesanderson::wes_palette("Darjeeling1", 2, "continuous")
valence = c("#119DA4", "#19647E")
plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

define functions

model_table = function(model) {
  model %>%
    broom.mixed::tidy(., conf.int = TRUE) %>%
    filter(effect == "fixed") %>%
    rename("SE" = std.error,
           "z" = statistic,
           "p" = p.value) %>%
    mutate(term = gsub("\\(Intercept\\)", "Intercept (self-oriented well-being)", term),
           term = gsub("log\\(RT\\)", "RT", term),
           term = gsub("constructill-being", "Construct (ill-being)", term),
           term = gsub("constructsocial well-being", "Construct (social well-being)", term),
           term = gsub(":", " x ", term),
           p = ifelse(p < .001, "< .001",
                      ifelse(p == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
           `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high),
           z = abs(round(z, 2))) %>%
    select(term, `b [95% CI]`, z, p) %>%
    kable() %>%
    kableExtra::kable_styling()
}

load data

task = read.csv("../../data/task_data.csv") %>%
  mutate(responseYN = as.factor(responseYN),
         construct = recode(construct, "well-being" = "self-oriented well-being",
                            "social" = "social well-being"),
         construct = factor(construct, levels = c("self-oriented well-being", "social well-being", "ill-being")),
         valence = factor(valence, levels = c("positive", "negative")))
betas = read.csv("../../data/neuro_data.csv")

exclude outliers and standardize

Subset self trials
Exclude outlier trials that are > 3 SDs from roi median across participants

trial_conditions = task %>%
  select(subjectID, trial, run, instruction)

betas_outlier = betas %>%
  left_join(., trial_conditions) %>%
  filter(instruction == "self") %>%
  group_by(roi) %>%
  mutate(median = median(meanPE, na.rm = TRUE),
         sd3 = 3*sd(meanPE, na.rm = TRUE),
         outlier = ifelse(meanPE > median + sd3 | meanPE < median - sd3, "yes", "no"))

## Joining, by = c("subjectID", "trial", "run")

betas_outlier %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = round((n / total) * 100, 2)) %>%
  kable(format = "pandoc")

outlier	n	total	percent
no	11900	11988	99.27
yes	88	11988	0.73

betas_ex = betas_outlier %>%
  filter(outlier == "no") %>%
  select(-c(median, sd3, outlier))

merge data and exclude participants

Motion exclusions: FP091
Technical failure: FP080, FP082
Non-compliance: FP021, FP049, FP085, FP121
Standardize within participant and ROI

data_ex = task %>%
  left_join(., betas_ex) %>%
  filter(!subjectID %in% c("FP091", "FP080", "FP082", "FP021", "FP049", "FP085", "FP121")) %>%
  filter(!is.na(meanPE)) %>%
  filter(instruction == "self") %>%
  group_by(roi, subjectID) %>%
  mutate(std = meanPE / sd(meanPE, na.rm = TRUE)) %>%
  select(-sdPE)

## Joining, by = c("trial", "subjectID", "run", "instruction")

remove missing data for MLM analysis

data_complete_mod = data_ex %>%
  select(-meanPE) %>%
  spread(roi, std) %>%
  select(subjectID, trial, run, construct, item, valence, response, responseYN, RT_original, RT, pgACC, vmPFC, VS) %>%
  na.omit()

item responses

percent endorsements

percents = data_complete_mod %>%
  group_by(construct, item, response) %>%
  summarize(n = n()) %>%
  spread(response, n) %>%
  mutate(no = ifelse(is.na(no), 0, no),
         yes = ifelse(is.na(yes), 0, yes),
         total = no + yes) %>%
  gather(response, n, no, yes) %>%
  mutate(percent = round((n/total) * 100, 1),
         n = sprintf("%s / %s", n, total)) %>%
  filter(response == "yes")

## `summarise()` has grouped output by 'construct', 'item'. You can override using
## the `.groups` argument.

percents %>%
  select(-response, -total) %>%
  arrange(desc(percent))  %>%
  kable() %>%
  kableExtra::kable_styling()

construct	item	n	percent
self-oriented well-being	capable	5 / 5	100.0
social well-being	have support systems	5 / 5	100.0
social well-being	well-liked	5 / 5	100.0
ill-being	tired	5 / 5	100.0
self-oriented well-being	interested in life	97 / 98	99.0
social well-being	loved	102 / 103	99.0
social well-being	trusted	101 / 104	97.1
social well-being	am cared for	98 / 102	96.1
social well-being	comfortable with friends	96 / 104	92.3
self-oriented well-being	happy	94 / 102	92.2
self-oriented well-being	live a meaningful life	95 / 103	92.2
self-oriented well-being	achieve goals	95 / 104	91.3
social well-being	have warm relationships	95 / 104	91.3
self-oriented well-being	hopeful	88 / 97	90.7
self-oriented well-being	live a purposeful life	92 / 102	90.2
self-oriented well-being	engaged in tasks	87 / 97	89.7
social well-being	have good friends	90 / 102	88.2
social well-being	have trusting relationships	87 / 99	87.9
self-oriented well-being	joyful	88 / 103	85.4
self-oriented well-being	like myself	84 / 103	81.6
social well-being	feel encouraged	82 / 101	81.2
self-oriented well-being	positive	79 / 98	80.6
self-oriented well-being	full of energy	4 / 5	80.0
self-oriented well-being	optimistic	4 / 5	80.0
social well-being	helped by friends	76 / 96	79.2
self-oriented well-being	satisfied with life	79 / 101	78.2
self-oriented well-being	proud	78 / 100	78.0
ill-being	stressed	67 / 100	67.0
social well-being	belong in social group	62 / 95	65.3
self-oriented well-being	cheerful	3 / 5	60.0
ill-being	anxious	57 / 104	54.8
ill-being	nervous	53 / 103	51.5
ill-being	overwhelmed	49 / 103	47.6
social well-being	lonely	31 / 102	30.4
ill-being	fearful	26 / 98	26.5
ill-being	unhappy	1 / 4	25.0
social well-being	isolated	24 / 102	23.5
ill-being	feel unimportant	18 / 103	17.5
ill-being	feel rejected	17 / 104	16.3
ill-being	sad	14 / 99	14.1
ill-being	angry	10 / 98	10.2
ill-being	feel useless	9 / 99	9.1
ill-being	unable to cope	8 / 104	7.7
ill-being	feel worthless	5 / 99	5.1
social well-being	alienated	0 / 5	0.0
social well-being	neglected	0 / 5	0.0
ill-being	depressed	0 / 4	0.0
ill-being	feel invisible	0 / 5	0.0

percent endorsements for improperly randomized items

percents %>%
  filter(total < 10) %>%
  select(-response, -total) %>%
  arrange(desc(construct)) %>%
  kable() %>%
  kableExtra::kable_styling()

construct	item	n	percent
ill-being	depressed	0 / 4	0
ill-being	feel invisible	0 / 5	0
ill-being	tired	5 / 5	100
ill-being	unhappy	1 / 4	25
social well-being	alienated	0 / 5	0
social well-being	have support systems	5 / 5	100
social well-being	neglected	0 / 5	0
social well-being	well-liked	5 / 5	100
self-oriented well-being	capable	5 / 5	100
self-oriented well-being	cheerful	3 / 5	60
self-oriented well-being	full of energy	4 / 5	80
self-oriented well-being	optimistic	4 / 5	80

plot

data_complete_mod %>%
  group_by(construct, item, responseYN) %>%
  summarize(n = n()) %>%
  spread(responseYN, n) %>%
  mutate(no = ifelse(is.na(no), 0, no),
         yes = ifelse(is.na(yes), 0, yes),
         total = no + yes) %>%
  gather(responseYN, n, no, yes) %>%
  mutate(percent = round((n/total) * 100, 1),
         n = sprintf("%s / %s", n, total)) %>%
  filter(responseYN == "yes") %>%
  ggplot(aes(percent, fill = construct)) +
  geom_density(color = NA, alpha = .8) +
  scale_fill_manual(values = constructs) +
  labs(x = "\n percent endorsed") +
  plot_aes

## `summarise()` has grouped output by 'construct', 'item'. You can override using
## the `.groups` argument.

well-being facet correlations

data_corr = data_complete_mod %>%
  mutate(responseNum = recode(responseYN, "yes" = 1, "no" = 0),
         construct = recode(construct, "self-oriented well-being" = "self",
                            "social well-being" = "social")) %>%
  group_by(construct, subjectID) %>%
  summarize(sum = sum(responseNum, na.rm = TRUE),
            n = n()) %>%
  mutate(percent = (sum / n) * 100) %>%
  select(subjectID, construct, percent) %>%
  spread(construct, percent)

## `summarise()` has grouped output by 'construct'. You can override using the
## `.groups` argument.

# correlations
cor.ci(data_corr[,2:4])

## Call:corCi(x = x, keys = keys, n.iter = n.iter, p = p, overlap = overlap, 
##     poly = poly, method = method, plot = plot, minlength = minlength, 
##     n = n)
## 
##  Coefficients and bootstrapped confidence intervals 
##           self  socil ill-b
## self       1.00            
## social     0.57  1.00      
## ill-being -0.39 -0.42  1.00
## 
##  scale correlations and bootstrapped confidence intervals 
##             lower.emp lower.norm estimate upper.norm upper.emp p
## self-socil       0.42       0.42     0.57       0.71      0.71 0
## self-ill-b      -0.50      -0.53    -0.39      -0.24     -0.26 0
## socil-ill-b     -0.60      -0.59    -0.42      -0.22     -0.25 0

ROI descriptives

means, SDs, and correlations

# correlations
pgACC_vmPFC = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, pgACC, vmPFC) %>%
   rmcorr::rmcorr(subjectID, pgACC, vmPFC, .)

pgACC_VS = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, pgACC, VS) %>%
   rmcorr::rmcorr(subjectID, pgACC, VS, .)

vmPFC_VS = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, vmPFC, VS) %>%
   rmcorr::rmcorr(subjectID, vmPFC, VS, .)

# means and SDs
mean_table = data_complete_mod %>%
  gather(roi, value, pgACC, vmPFC, VS) %>%
  group_by(roi) %>%
  summarize(M = round(mean(value, na.rm = TRUE), 2),
            SD = round(sd(value, na.rm = TRUE), 2)) %>%
  mutate(order = ifelse(roi == "pgACC", 1,
                 ifelse(roi == "vmPFC", 2, 3)),
         roi = ifelse(roi == "pgACC", "1. pgACC",
               ifelse(roi == "vmPFC", "2. vmPFC", "3. VS"))) %>%
  rename("ROI" = roi) %>%
  arrange(order) %>%
  select(-order)

# table
corr_table = data.frame(ROI = c("1. pgACC", "2. vmPFC", "3. VS"),
           `1` = c("--", paste0(round(pgACC_vmPFC[[1]],2), " [",round(pgACC_vmPFC[[4]][1],2),", ",round(pgACC_vmPFC[[4]][2],2),"]"),
                   paste0(round(pgACC_VS[[1]],2), " [",round(pgACC_VS[[4]][1],2),", ",round(pgACC_VS[[4]][2],2),"]")),
           `2` = c("", "--", paste0(round(vmPFC_VS[[1]],2), " [",round(vmPFC_VS[[4]][1],2),", ",round(vmPFC_VS[[4]][2],2),"]")),
           `3` = c("", "", "--"),
           check.names = FALSE)

mean_table %>%
  full_join(., corr_table) %>%
  knitr::kable(format = "pandoc", caption = "Repeated Measures Correlations Between ROIs (n = 3694). All correlations are statistically significant, p < .001. 95% confidence intervals are bracketed")

## Joining, by = "ROI"

Repeated Measures Correlations Between ROIs (n = 3694). All correlations are statistically significant, p < .001. 95% confidence intervals are bracketed
ROI	M	SD	1	2	3
1. pgACC	0.04	1.03	–
2. vmPFC	0.15	1.02	0.84 [0.83, 0.85]	–
3. VS	0.21	1.03	0.42 [0.39, 0.45]	0.45 [0.43, 0.48]	–

RT descriptives

across facets and responses

data_complete_mod %>%
  ungroup() %>%
  summarize(`mean RT` = mean(RT_original, na.rm = TRUE),
            `median RT` = median(RT_original, na.rm = TRUE),
            SD = sd(RT_original, na.rm = TRUE)) %>%
    kable(digits = 2) %>%
    kableExtra::kable_styling()

mean RT	median RT	SD
1.62	1.48	0.61

across responses

data_complete_mod %>%
  group_by(construct) %>%
  summarize(`mean RT` = mean(RT_original, na.rm = TRUE),
            `median RT` = median(RT_original, na.rm = TRUE),
            SD = sd(RT_original, na.rm = TRUE)) %>%
    kable(digits = 2) %>%
    kableExtra::kable_styling()

construct	mean RT	median RT	SD
self-oriented well-being	1.53	1.40	0.57
social well-being	1.61	1.44	0.62
ill-being	1.72	1.58	0.62

data_complete_mod %>%
  ggplot(aes(RT_original, fill = construct)) +
  geom_density(color = NA, alpha = .75) +
  scale_fill_manual(values = constructs) +
  labs(x = "\nRT (s)") +
  plot_aes

by response

data_complete_mod %>%
  group_by(construct, responseYN) %>%
  summarize(`mean RT` = mean(RT_original, na.rm = TRUE)) %>%
  spread(responseYN, `mean RT`) %>%
    kable() %>%
    kableExtra::kable_styling()

## `summarise()` has grouped output by 'construct'. You can override using the
## `.groups` argument.

construct	no	yes
self-oriented well-being	2.015276	1.455663
social well-being	2.113079	1.530455
ill-being	1.709059	1.746566

model testing diffferences

model_rt = lmer(RT_original ~ responseYN * construct + (1 | subjectID),
                data = data_complete_mod)

summary(model_rt)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT_original ~ responseYN * construct + (1 | subjectID)
##    Data: data_complete_mod
## 
## REML criterion at convergence: 5910.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8363 -0.6517 -0.1830  0.4757  5.8536 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.07652  0.2766  
##  Residual              0.26907  0.5187  
## Number of obs: 3694, groups:  subjectID, 105
## 
## Fixed effects:
##                                            Estimate Std. Error         df
## (Intercept)                                 1.93039    0.05072  894.26616
## responseYNyes                              -0.46016    0.04614 3632.65369
## constructsocial well-being                  0.08782    0.05777 3599.67437
## constructill-being                         -0.18627    0.04660 3629.89748
## responseYNyes:constructsocial well-being   -0.00823    0.06222 3602.24250
## responseYNyes:constructill-being            0.38339    0.05826 3646.40501
##                                          t value             Pr(>|t|)    
## (Intercept)                               38.062 < 0.0000000000000002 ***
## responseYNyes                             -9.973 < 0.0000000000000002 ***
## constructsocial well-being                 1.520                0.129    
## constructill-being                        -3.997      0.0000653894996 ***
## responseYNyes:constructsocial well-being  -0.132                0.895    
## responseYNyes:constructill-being           6.581      0.0000000000535 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) rspnYN cnwll- cnstr- rYN:w-
## responsYNys -0.794                            
## cnstrctwll- -0.607  0.670                     
## cnstrctll-b -0.784  0.870  0.661              
## rspnsYNy:w-  0.566 -0.712 -0.932 -0.616       
## rspnsYNys:-  0.643 -0.810 -0.532 -0.823  0.565

predicted_sub = ggpredict(model_rt, terms = c("responseYN", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_rt, terms = c("responseYN", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  #geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .25) +
  geom_line(aes(color = group, group = group), size = 1.5) +
  geom_pointrange(aes(color = group, ymin = conf.low, ymax = conf.high)) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nresponse", y = "reaction time (s)\n") + 
  plot_aes

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

RT models

These are the primary analyses reported in the main text of the manuscript, controlling for reaction time.

run models

compare models

Determine the best fitting model.

model_0_rt = Base model including including no ROIs

model_1_rt = Main effects model including ROIs

model_2_rt = Interaction model including ROIs and their interactions with well-being construct

model_3_rt = Interaction model including ROIs and their interactions with valence rather than well-being construct

model_0_rt = glmer(responseYN ~ construct * RT + 
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_1_rt = glmer(responseYN ~ construct * RT + pgACC + vmPFC + VS + 
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_2_rt = glmer(responseYN ~ RT*construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_3_rt = glmer(responseYN ~ RT*valence + pgACC*valence + vmPFC*valence + VS*valence +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

# compare models
anova(model_0_rt, model_1_rt, model_2_rt, model_3_rt)

calculate R2

Compare the model R2 for the best fitting compared to the base model.

performance::r2(model_0_rt)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.512
##      Marginal R2: 0.479

performance::r2(model_2_rt)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.532
##      Marginal R2: 0.498

check variance inflation in the best fitting model

# vif
car::vif(model_2_rt)

##                      GVIF Df GVIF^(1/(2*Df))
## RT               3.985384  1        1.996343
## construct        1.514415  2        1.109331
## pgACC           13.575127  1        3.684444
## vmPFC           13.404569  1        3.661225
## VS               4.932306  1        2.220880
## RT:construct     4.236780  2        1.434693
## construct:pgACC 40.780615  2        2.527048
## construct:vmPFC 48.147820  2        2.634172
## construct:VS     6.643231  2        1.605443

model summary

summary(model_2_rt)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ RT * construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3082.3   3181.7  -1525.1   3050.3     3678 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1633 -0.5202  0.2040  0.3958  3.1031 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.2382   0.488   
## Number of obs: 3694, groups:  subjectID, 105
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       1.89473    0.11189  16.933
## RT                               -3.15503    0.29977 -10.525
## constructsocial well-being        0.08607    0.14486   0.594
## constructill-being               -2.95012    0.12509 -23.584
## pgACC                             0.44782    0.18871   2.373
## vmPFC                            -0.39372    0.18837  -2.090
## VS                                0.10966    0.11200   0.979
## RT:constructsocial well-being    -0.03285    0.39988  -0.082
## RT:constructill-being             3.16002    0.35531   8.894
## constructsocial well-being:pgACC  0.06793    0.25779   0.264
## constructill-being:pgACC         -0.38844    0.22428  -1.732
## constructsocial well-being:vmPFC -0.23103    0.26387  -0.876
## constructill-being:vmPFC          0.75862    0.22610   3.355
## constructsocial well-being:VS    -0.03716    0.15559  -0.239
## constructill-being:VS            -0.18072    0.13704  -1.319
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## RT                               < 0.0000000000000002 ***
## constructsocial well-being                   0.552418    
## constructill-being               < 0.0000000000000002 ***
## pgACC                                        0.017643 *  
## vmPFC                                        0.036604 *  
## VS                                           0.327535    
## RT:constructsocial well-being                0.934521    
## RT:constructill-being            < 0.0000000000000002 ***
## constructsocial well-being:pgACC             0.792143    
## constructill-being:pgACC                     0.083287 .  
## constructsocial well-being:vmPFC             0.381276    
## constructill-being:vmPFC                     0.000793 ***
## constructsocial well-being:VS                0.811232    
## constructill-being:VS                        0.187246    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

model table

model_table(model_2_rt)

term	b [95% CI]	z	p
Intercept (self-oriented well-being)	1.89 [1.68, 2.11]	16.93	< .001
RT	-3.16 [-3.74, -2.57]	10.52	< .001
Construct (social well-being)	0.09 [-0.20, 0.37]	0.59	.552
Construct (ill-being)	-2.95 [-3.20, -2.70]	23.58	< .001
pgACC	0.45 [0.08, 0.82]	2.37	.018
vmPFC	-0.39 [-0.76, -0.02]	2.09	.037
VS	0.11 [-0.11, 0.33]	0.98	.328
RT x Construct (social well-being)	-0.03 [-0.82, 0.75]	0.08	.935
RT x Construct (ill-being)	3.16 [2.46, 3.86]	8.89	< .001
Construct (social well-being) x pgACC	0.07 [-0.44, 0.57]	0.26	.792
Construct (ill-being) x pgACC	-0.39 [-0.83, 0.05]	1.73	.083
Construct (social well-being) x vmPFC	-0.23 [-0.75, 0.29]	0.88	.381
Construct (ill-being) x vmPFC	0.76 [0.32, 1.20]	3.36	< .001
Construct (social well-being) x VS	-0.04 [-0.34, 0.27]	0.24	.811
Construct (ill-being) x VS	-0.18 [-0.45, 0.09]	1.32	.187

marginal probabilities

all

ggeffect(model_2_rt)

## $RT
## # Predicted probabilities of responseYN
## 
##    RT | Predicted |       95% CI
## --------------------------------
## -1.50 |      0.98 | [0.97, 0.99]
## -1.00 |      0.95 | [0.94, 0.97]
## -0.50 |      0.88 | [0.85, 0.90]
##  0.00 |      0.72 | [0.69, 0.74]
##  0.50 |      0.47 | [0.42, 0.52]
##  1.00 |      0.23 | [0.18, 0.30]
## 
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.90 | [0.88, 0.92]
## social well-being        |      0.91 | [0.88, 0.93]
## ill-being                |      0.27 | [0.24, 0.30]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.29 | [0.11, 0.56]
##    -4 |      0.45 | [0.27, 0.64]
##    -2 |      0.62 | [0.52, 0.70]
##     0 |      0.76 | [0.73, 0.79]
##     2 |      0.86 | [0.80, 0.90]
##     4 |      0.93 | [0.85, 0.96]
##     6 |      0.96 | [0.88, 0.99]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -4 |      0.89 | [0.77, 0.95]
##    -2 |      0.84 | [0.76, 0.89]
##     0 |      0.77 | [0.74, 0.79]
##     2 |      0.68 | [0.60, 0.76]
##     4 |      0.58 | [0.39, 0.75]
##     6 |      0.47 | [0.22, 0.74]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -4 |      0.73 | [0.62, 0.82]
## -2 |      0.75 | [0.69, 0.80]
##  0 |      0.76 | [0.73, 0.79]
##  2 |      0.77 | [0.72, 0.82]
##  4 |      0.79 | [0.70, 0.86]
##  6 |      0.80 | [0.66, 0.89]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2_rt"

pgACC

ggeffect(model_2_rt, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.9013345	0.1186334	0.8786410	0.9201699	0	0.90 [0.88, 0.92]
social well-being	0.9055447	0.1222800	0.8829581	0.9241470	0	0.91 [0.88, 0.92]
ill-being	0.2655698	0.0838974	0.2347550	0.2988497	0	0.27 [0.23, 0.30]
self-oriented well-being	0.9346223	0.2305713	0.9009700	0.9573798	1	0.93 [0.90, 0.96]
social well-being	0.9413741	0.2312392	0.9107595	0.9619252	1	0.94 [0.91, 0.96]
ill-being	0.2773118	0.1447915	0.2241538	0.3375907	1	0.28 [0.22, 0.34]

vmPFC

ggeffect(model_2_rt, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.9078067	0.1249982	0.8851522	0.9263642	0	0.91 [0.89, 0.93]
social well-being	0.9146557	0.1344663	0.8917078	0.9331056	0	0.91 [0.89, 0.93]
ill-being	0.2557264	0.0861418	0.2249350	0.2891602	0	0.26 [0.22, 0.29]
self-oriented well-being	0.8691456	0.1913317	0.8203060	0.9062283	1	0.87 [0.82, 0.91]
social well-being	0.8515878	0.1704104	0.8042561	0.8890533	1	0.85 [0.80, 0.89]
ill-being	0.3310591	0.1367623	0.2745920	0.3928496	1	0.33 [0.27, 0.39]

VS

ggeffect(model_2_rt, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.9008684	0.1201313	0.8777676	0.9200013	0	0.90 [0.88, 0.92]
social well-being	0.9059780	0.1228265	0.8833714	0.9245769	0	0.91 [0.88, 0.92]
ill-being	0.2688924	0.0857897	0.2371454	0.3032002	0	0.27 [0.24, 0.30]
self-oriented well-being	0.9102403	0.1535243	0.8824326	0.9319778	1	0.91 [0.88, 0.93]
social well-being	0.9119746	0.1598452	0.8833657	0.9340897	1	0.91 [0.88, 0.93]
ill-being	0.2551545	0.1038266	0.2184360	0.2957100	1	0.26 [0.22, 0.30]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2_rt, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2_rt, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2_rt, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2_rt, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2_rt, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2_rt, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

RT

vals = seq(round(min(data_complete_mod$RT),1), round(max(data_complete_mod$RT), 1), .1)

predicted_sub = ggpredict(model_2_rt, terms = c("RT [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_2_rt, terms = c("RT [vals]", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 1.5) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nlog-transformed grand-mean centered RT (s)", y = "probability of responding yes\n") + 
  plot_aes

no RT model

These are the original preregistered analyses not controlling for reaction time, reported in supplementary material.

run model

model_2 = glmer(responseYN ~ construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model summary

summary(model_2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3355.6   3436.4  -1664.8   3329.6     3681 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6442 -0.5280  0.3366  0.4098  3.0378 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.1913   0.4373  
## Number of obs: 3694, groups:  subjectID, 105
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       2.01006    0.10105  19.893
## constructsocial well-being       -0.09374    0.12688  -0.739
## constructill-being               -3.05942    0.11701 -26.146
## pgACC                             0.24530    0.15334   1.600
## vmPFC                            -0.22700    0.15458  -1.469
## VS                                0.12710    0.09061   1.403
## constructsocial well-being:pgACC  0.01747    0.21058   0.083
## constructill-being:pgACC         -0.19387    0.19550  -0.992
## constructsocial well-being:vmPFC -0.07919    0.21652  -0.366
## constructill-being:vmPFC          0.59383    0.19900   2.984
## constructsocial well-being:VS    -0.07809    0.12521  -0.624
## constructill-being:VS            -0.19412    0.12025  -1.614
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## constructsocial well-being                    0.46005    
## constructill-being               < 0.0000000000000002 ***
## pgACC                                         0.10966    
## vmPFC                                         0.14196    
## VS                                            0.16071    
## constructsocial well-being:pgACC              0.93388    
## constructill-being:pgACC                      0.32137    
## constructsocial well-being:vmPFC              0.71457    
## constructill-being:vmPFC                      0.00285 ** 
## constructsocial well-being:VS                 0.53281    
## constructill-being:VS                         0.10645    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnwll- cnstr- pgACC  vmPFC  VS     cw-:AC c-:ACC cw-:PF
## cnstrctwll- -0.627                                                        
## cnstrctll-b -0.725  0.543                                                 
## pgACC        0.210 -0.164 -0.182                                          
## vmPFC       -0.179  0.139  0.155 -0.804                                   
## VS          -0.079  0.060  0.068 -0.087 -0.197                            
## cwll-bn:ACC -0.156  0.261  0.137 -0.718  0.577  0.062                     
## cnstrc-:ACC -0.161  0.129  0.184 -0.773  0.621  0.067  0.562              
## cwll-bn:PFC  0.128 -0.286 -0.112  0.565 -0.703  0.139 -0.799 -0.443       
## cnstrc-:PFC  0.142 -0.109 -0.188  0.615 -0.766  0.151 -0.448 -0.793  0.546
## cwll-bng:VS  0.053 -0.057 -0.045  0.060  0.140 -0.712 -0.075 -0.049 -0.215
## cnstrct-:VS  0.058 -0.044 -0.128  0.066  0.145 -0.745 -0.048 -0.078 -0.104
##             c-:PFC cw-:VS
## cnstrctwll-              
## cnstrctll-b              
## pgACC                    
## vmPFC                    
## VS                       
## cwll-bn:ACC              
## cnstrc-:ACC              
## cwll-bn:PFC              
## cnstrc-:PFC              
## cwll-bng:VS -0.109       
## cnstrct-:VS -0.208  0.537

model table

model_table(model_2)

term	b [95% CI]	z	p
Intercept (self-oriented well-being)	2.01 [1.81, 2.21]	19.89	< .001
Construct (social well-being)	-0.09 [-0.34, 0.15]	0.74	.460
Construct (ill-being)	-3.06 [-3.29, -2.83]	26.15	< .001
pgACC	0.25 [-0.06, 0.55]	1.60	.110
vmPFC	-0.23 [-0.53, 0.08]	1.47	.142
VS	0.13 [-0.05, 0.30]	1.40	.161
Construct (social well-being) x pgACC	0.02 [-0.40, 0.43]	0.08	.934
Construct (ill-being) x pgACC	-0.19 [-0.58, 0.19]	0.99	.321
Construct (social well-being) x vmPFC	-0.08 [-0.50, 0.35]	0.37	.715
Construct (ill-being) x vmPFC	0.59 [0.20, 0.98]	2.98	.003
Construct (social well-being) x VS	-0.08 [-0.32, 0.17]	0.62	.533
Construct (ill-being) x VS	-0.19 [-0.43, 0.04]	1.61	.106

marginal probabilities

all

ggeffect(model_2)

## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.88 | [0.86, 0.90]
## social well-being        |      0.87 | [0.85, 0.89]
## ill-being                |      0.27 | [0.24, 0.30]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.46 | [0.24, 0.70]
##    -4 |      0.55 | [0.39, 0.71]
##    -2 |      0.64 | [0.56, 0.72]
##     0 |      0.72 | [0.70, 0.75]
##     2 |      0.79 | [0.73, 0.84]
##     4 |      0.85 | [0.74, 0.91]
##     6 |      0.89 | [0.75, 0.96]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -4 |      0.77 | [0.62, 0.87]
##    -2 |      0.75 | [0.67, 0.81]
##     0 |      0.73 | [0.70, 0.75]
##     2 |      0.70 | [0.63, 0.77]
##     4 |      0.68 | [0.53, 0.80]
##     6 |      0.65 | [0.42, 0.84]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -4 |      0.69 | [0.59, 0.78]
## -2 |      0.71 | [0.65, 0.76]
##  0 |      0.72 | [0.70, 0.75]
##  2 |      0.74 | [0.69, 0.78]
##  4 |      0.75 | [0.67, 0.82]
##  6 |      0.76 | [0.64, 0.85]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2"

pgACC

ggeffect(model_2, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.8811189	0.0989369	0.8592598	0.8999781	0	0.88 [0.86, 0.90]
social well-being	0.8677987	0.0961826	0.8446340	0.8879677	0	0.87 [0.84, 0.89]
ill-being	0.2670746	0.0795909	0.2376673	0.2986948	0	0.27 [0.24, 0.30]
self-oriented well-being	0.9045098	0.1837014	0.8685640	0.9314013	1	0.90 [0.87, 0.93]
social well-being	0.8951448	0.1805836	0.8569846	0.9240255	1	0.90 [0.86, 0.92]
ill-being	0.2772620	0.1410630	0.2253839	0.3359031	1	0.28 [0.23, 0.34]

vmPFC

ggeffect(model_2, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.8855238	0.1026127	0.8635030	0.9043853	0	0.89 [0.86, 0.90]
social well-being	0.8739828	0.1036531	0.8498577	0.8947116	0	0.87 [0.85, 0.89]
ill-being	0.2570822	0.0816078	0.2277362	0.2887955	0	0.26 [0.23, 0.29]
self-oriented well-being	0.8604231	0.1619427	0.8177841	0.8943731	1	0.86 [0.82, 0.89]
social well-being	0.8362314	0.1466111	0.7929990	0.8718932	1	0.84 [0.79, 0.87]
ill-being	0.3330630	0.1345046	0.2772812	0.3939497	1	0.33 [0.28, 0.39]

VS

ggeffect(model_2, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.8793589	0.0998866	0.8570007	0.8986349	0	0.88 [0.86, 0.90]
social well-being	0.8678097	0.0966799	0.8445186	0.8880742	0	0.87 [0.84, 0.89]
ill-being	0.2701836	0.0813785	0.2399069	0.3027594	0	0.27 [0.24, 0.30]
self-oriented well-being	0.8922061	0.1258366	0.8660923	0.9137346	1	0.89 [0.87, 0.91]
social well-being	0.8733311	0.1246071	0.8437650	0.8979786	1	0.87 [0.84, 0.90]
ill-being	0.2571747	0.1004826	0.2213791	0.2965537	1	0.26 [0.22, 0.30]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

item random effect model

This is a model fit in response to reviewers, included in supplementary material.

run model

model_2_rt_item = glmer(responseYN ~ RT*construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID) + (1 | item), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model summary

summary(model_2_rt_item)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ RT * construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID) + (1 | item)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   2786.0   2891.6  -1376.0   2752.0     3677 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.9511  -0.3244   0.1677   0.3614   4.6525 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.3501   0.5917  
##  item      (Intercept) 1.0609   1.0300  
## Number of obs: 3694, groups:  subjectID, 105; item, 48
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       1.96805    0.30826   6.384
## RT                               -3.51299    0.31973 -10.988
## constructsocial well-being        0.41306    0.43415   0.951
## constructill-being               -3.24427    0.41970  -7.730
## pgACC                             0.39922    0.19094   2.091
## vmPFC                            -0.36949    0.19099  -1.935
## VS                                0.13200    0.11471   1.151
## RT:constructsocial well-being     0.34789    0.42647   0.816
## RT:constructill-being             4.05498    0.39109  10.368
## constructsocial well-being:pgACC  0.13295    0.26468   0.502
## constructill-being:pgACC         -0.41633    0.23229  -1.792
## constructsocial well-being:vmPFC -0.22625    0.27384  -0.826
## constructill-being:vmPFC          0.66764    0.23558   2.834
## constructsocial well-being:VS    -0.07297    0.16141  -0.452
## constructill-being:VS            -0.12616    0.14516  -0.869
##                                              Pr(>|z|)    
## (Intercept)                        0.0000000001721625 ***
## RT                               < 0.0000000000000002 ***
## constructsocial well-being                     0.3414    
## constructill-being                 0.0000000000000107 ***
## pgACC                                          0.0365 *  
## vmPFC                                          0.0530 .  
## VS                                             0.2498    
## RT:constructsocial well-being                  0.4146    
## RT:constructill-being            < 0.0000000000000002 ***
## constructsocial well-being:pgACC               0.6155    
## constructill-being:pgACC                       0.0731 .  
## constructsocial well-being:vmPFC               0.4087    
## constructill-being:vmPFC                       0.0046 ** 
## constructsocial well-being:VS                  0.6512    
## constructill-being:VS                          0.3848    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

model table

model_table(model_2_rt_item)

term	b [95% CI]	z	p
Intercept (self-oriented well-being)	1.97 [1.36, 2.57]	6.38	< .001
RT	-3.51 [-4.14, -2.89]	10.99	< .001
Construct (social well-being)	0.41 [-0.44, 1.26]	0.95	.341
Construct (ill-being)	-3.24 [-4.07, -2.42]	7.73	< .001
pgACC	0.40 [0.02, 0.77]	2.09	.037
vmPFC	-0.37 [-0.74, 0.00]	1.93	.053
VS	0.13 [-0.09, 0.36]	1.15	.250
RT x Construct (social well-being)	0.35 [-0.49, 1.18]	0.82	.415
RT x Construct (ill-being)	4.05 [3.29, 4.82]	10.37	< .001
Construct (social well-being) x pgACC	0.13 [-0.39, 0.65]	0.50	.615
Construct (ill-being) x pgACC	-0.42 [-0.87, 0.04]	1.79	.073
Construct (social well-being) x vmPFC	-0.23 [-0.76, 0.31]	0.83	.409
Construct (ill-being) x vmPFC	0.67 [0.21, 1.13]	2.83	.005
Construct (social well-being) x VS	-0.07 [-0.39, 0.24]	0.45	.651
Construct (ill-being) x VS	-0.13 [-0.41, 0.16]	0.87	.385

marginal probabilities

all

ggeffect(model_2_rt_item)

## $RT
## # Predicted probabilities of responseYN
## 
##    RT | Predicted |       95% CI
## --------------------------------
## -1.50 |      0.98 | [0.97, 0.99]
## -1.00 |      0.96 | [0.93, 0.97]
## -0.50 |      0.88 | [0.84, 0.92]
##  0.00 |      0.73 | [0.66, 0.80]
##  0.50 |      0.50 | [0.40, 0.60]
##  1.00 |      0.26 | [0.18, 0.37]
## 
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.91 | [0.85, 0.95]
## social well-being        |      0.94 | [0.89, 0.96]
## ill-being                |      0.22 | [0.13, 0.33]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.36 | [0.14, 0.66]
##    -4 |      0.50 | [0.30, 0.71]
##    -2 |      0.65 | [0.52, 0.76]
##     0 |      0.77 | [0.70, 0.83]
##     2 |      0.86 | [0.79, 0.92]
##     4 |      0.92 | [0.83, 0.97]
##     6 |      0.96 | [0.86, 0.99]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -4 |      0.90 | [0.77, 0.96]
##    -2 |      0.85 | [0.76, 0.91]
##     0 |      0.78 | [0.71, 0.84]
##     2 |      0.70 | [0.58, 0.79]
##     4 |      0.60 | [0.38, 0.78]
##     6 |      0.49 | [0.21, 0.77]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -4 |      0.73 | [0.58, 0.83]
## -2 |      0.75 | [0.66, 0.83]
##  0 |      0.77 | [0.70, 0.83]
##  2 |      0.80 | [0.72, 0.86]
##  4 |      0.82 | [0.71, 0.89]
##  6 |      0.84 | [0.69, 0.92]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2_rt_item"

pgACC

ggeffect(model_2_rt_item, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.9116221	0.3103964	0.8488022	0.9498825	0	0.91 [0.85, 0.95]
social well-being	0.9346029	0.3199374	0.8841754	0.9639701	0	0.93 [0.88, 0.96]
ill-being	0.2155277	0.2958810	0.1333292	0.3291558	0	0.22 [0.13, 0.33]
self-oriented well-being	0.9389360	0.3698387	0.8816297	0.9694598	1	0.94 [0.88, 0.97]
social well-being	0.9605250	0.3804260	0.9202824	0.9808747	1	0.96 [0.92, 0.98]
ill-being	0.2126499	0.3240925	0.1251833	0.3376428	1	0.21 [0.13, 0.34]

vmPFC

ggeffect(model_2_rt_item, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.9170690	0.3129186	0.8569107	0.9533132	0	0.92 [0.86, 0.95]
social well-being	0.9408932	0.3251703	0.8938003	0.9678544	0	0.94 [0.89, 0.97]
ill-being	0.2081341	0.2967131	0.1281119	0.3198061	0	0.21 [0.13, 0.32]
self-oriented well-being	0.8842881	0.3460413	0.7950193	0.9377256	1	0.88 [0.80, 0.94]
social well-being	0.8976827	0.3455053	0.8167617	0.9452626	1	0.90 [0.82, 0.95]
ill-being	0.2615243	0.3191904	0.1592709	0.3983226	1	0.26 [0.16, 0.40]

VS

ggeffect(model_2_rt_item, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()

x	predicted	std.error	conf.low	conf.high	group	b [95% CI]
self-oriented well-being	0.9106727	0.3113140	0.8470583	0.9494078	0	0.91 [0.85, 0.95]
social well-being	0.9351202	0.3203609	0.8849579	0.9642925	0	0.94 [0.88, 0.96]
ill-being	0.2152117	0.2968624	0.1328915	0.3291677	0	0.22 [0.13, 0.33]
self-oriented well-being	0.9208444	0.3249766	0.8601977	0.9565124	1	0.92 [0.86, 0.96]
social well-being	0.9386109	0.3371729	0.8875865	0.9673279	1	0.94 [0.89, 0.97]
ill-being	0.2161998	0.3029340	0.1321948	0.3330969	1	0.22 [0.13, 0.33]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2_rt_item, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2_rt_item, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2_rt_item, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2_rt_item, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2_rt_item, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2_rt_item, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

RT

vals = seq(round(min(data_complete_mod$RT),1), round(max(data_complete_mod$RT), 1), .1)

predicted_sub = ggpredict(model_2_rt_item, terms = c("RT [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_2_rt_item, terms = c("RT [vals]", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 1.5) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nlog-transformed grand-mean centered RT (s)", y = "probability of responding yes\n") + 
  plot_aes

Trial-level analyses

Dani Cosme

2023-07-18

load packages

define aesthetics

define functions

load data

exclude outliers and standardize

merge data and exclude participants

remove missing data for MLM analysis

item responses

percent endorsements

percent endorsements for improperly randomized items

plot

well-being facet correlations

ROI descriptives

means, SDs, and correlations

RT descriptives

across facets and responses

across responses

by response

model testing diffferences

RT models

run models

compare models

calculate R2

check variance inflation in the best fitting model

model summary

model table

marginal probabilities

all

pgACC

vmPFC

VS

plot predicted effects

by construct

by ROI

RT

no RT model

run model

model summary

model table

marginal probabilities

all

pgACC

vmPFC

VS

plot predicted effects

by construct

by ROI

item random effect model

run model

model summary

model table

marginal probabilities

all

pgACC

vmPFC

VS

plot predicted effects

by construct

by ROI

RT